include("dg.jl")
using Plots

function main()
    # Order of method (m), number of elements (N)
    m = 4; N = 8

    # Set problem parameters
    FinalTime = 1.5; CFL = 0.25

    # Define spatial grid
    xL = -1.0; xR = 1.0; L = xR - xL; h = L/N

    # Generate mesh
    r = legendregausslobatto(m); VX = [xL + i * h for i = 0:N]
    x = ones(m + 1, 1) * VX[1:N]' + 0.5 * (r .+ 1) * (VX[2:N+1] - VX[1:N])'

    vmapM, vmapP, vmapI, vmapO = geomatry(m,N)

    D, rx, LIFT, Fscale = matrics(m,r,x)

    # Define initial conditions
    ϵ = 0.1; u = -tanh.((x .+ 0.5)/(2*ϵ)) .+ 1.0

    # Compute smallest spatial scale timestep
    time = 0; resu = zeros(m+1, N); xmin = minimum(abs.(x[1,:] - x[2,:]));
    umax = maximum(maximum(abs.(u), dims=2))
    dt = CFL* min(xmin/umax,xmin^2/sqrt(ϵ)); Nsteps = ceil(FinalTime/dt); dt = FinalTime/Nsteps

    # define rk coefficients
    rk4a = [0.0,
            -567301805773.0/1357537059087.0,
            -2404267990393.0/2016746695238.0,
            -3550918686646.0/2091501179385.0,
            -1275806237668.0/842570457699.0];
    rk4b = [1432997174477.0/9575080441755.0,
            5161836677717.0/13612068292357.0,
            1720146321549.0/2090206949498.0,
            3134564353537.0/4481467310338.0,
            2277821191437.0/14882151754819.0]
    rk4c = [0.0,
            1432997174477.0/9575080441755.0,
            2526269341429.0/6820363962896.0,
            2006345519317.0/3224310063776.0,
            2802321613138.0/2924317926251.0];

    # time step loop
    for tstep = 1:Nsteps
        for i  = 1:5
            time_local = time + rk4c[i]*dt
            rhsu = rhs(time_local,ϵ,xL,xR,u,h,m,N,D,vmapM,vmapP,vmapI,vmapO,rx, LIFT, Fscale)
            resu = rk4a[i]*resu + dt*rhsu
            u = u + rk4b[i]*resu
        end
        time = time+dt
    end
    x, u
end
